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A SIMPLE DISCRETIZATION SCHEME 



Abstract. A discretization scheme for nonnegative diffusion processes is proposed 
and the convergence of the corresponding sequence of approximate processes is proved 
using the martingale problem framework. Motivations for this scheme come typically 
from finance, especially for path-dependent option pricing. The scheme is simple: 
one only needs to find a nonnegative distribution whose mean and variance satisfy 
a simple condition to apply it. Then, for virtually any (path-dependent) payoff, 
Monte Carlo option prices obtained from this scheme will converge to the theoretical 
price. Examples of models and diffusion processes for which the scheme applies are 
provided. 



AMS Subject Classifications: 60J35, 65C30, 60H35, 91B24. 



Keywords: Euler discretization schemes, nonnegativity preservation, diffusion processes, Markov 
chains, martingale problem, convergence in distribution, interest rate models, stochastic volatil- 
ity models, path-dependent options. 



1. Introduction 



The Cox-Ingersoll-Ross (CIR) process, also known as the mean-reverting square-root diffusion. 



was introduced by Cox et al. (1985) for interest rates modeling. It now has other financial ap- 



plications, for example in Heston's stochastic volatility model (Heston, 1993), where it plays the 



role of the squared volatility. This process is the solution to the following stochastic differential 
equation: 



(1) 



dX{t) = K{f3 - X{t)) dt + i^^/X{t) dW{t), X{0) = xo, 



where {W{t))t>o is a one-dimensional standard Brownian motion, k, (3 and v are strictly positive 
constants, and the initial value satisfies xq > 0. It is known that this process stays nonnegative. 

If one wants to construct a discrete-time approximation of {X{t))t>Q defined by ([T]), a standard 
Euler-Maruyama scheme cannot be applied directly. Indeed, for a time-step of size 1/n, the 
approximating process {Yn{k))k>Q would be given by 1^(0) = xq, and 



(2) 



Yn{k + l) 



n 



Yn{k)) + 



n 



VYn(k)Zt\ 



(n) 

for A; = 0, 1, 2, ... , where the s are random variables with the standard normal distribution. 
If Yn(k) is nonnegative for a given k, then ^ correctly defines Yn{k + 1), whose value then 
has a non-zero probability of being negative, li Yn{k + 1) actually happens to be negative, the 
next iteration cannot be achieved since it would involve taking its square root. Thus, one can 
iterate (12]) only until the first k for which Yn{k) < 0. Several ways have been proposed for 

(n) 

avoiding this problem. For example, one may simulate values of until the right-hand side 
of ([2| is nonnegative, and then set y„(A; -|- 1) to be this value, but this results in a scheme for 
which the number of steps needed to generate a sample of a given size is random. One could 
also use either of these schemes instead of ([2|: 



(bl) Ynik + 1) 
(b2) Yn{k+1) 



Yn{k) + ^i(3-Y4k)) + 

Yn{k) + l^{P-{Yn{k)) + ) 



u 



in). 



in). 
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(b3) Yn{k + 1) =1 Yn{k) + ^(13 - Yn{k)) + 4"^ 
(b4) Yn{k + 1) = Yn{k) + - Yn{k)) + , 

where x'^ := max(x,0), x G M. All these schemes are well defined, but (bl), (b2) and (b4) 
generate processes {Yn{k))k>o whose values are not necessarily nonnegative. In quantitative 
finance, this is often perceived as a drawback when approximating positive quantities, such as 
interest rates, stock prices and volatilities. For more details on Euler-Maruyama discretization 
schemes for diffusions, the reader is referred to Glasserman (2004) and Kloeden and Platen 
( 1992] ). 

A natural yet crucial question when using a discretization scheme is: does it converge to the 
process we wish to approximate as the time-step decreases to zero? Classical theory mostly 
deals with diffusions with Lipschitz coefficients, excluding the CIR process (whose diffusion 
coefficient is not Lipschitz). [Deelstra and Delbaen (1998) establish the strong convergence of 
the scheme (bl), in a framework where the mean reversion parameter /3 may be a stochastic 
process. More recently. Bossy and Diop (2007) and Berkaoui et al. (2008) studied the weak 
and strong convergence of the scheme (b3), under the more general setting where the diffusion 
coefficient in ([T| is replaced with j>(X(t))", for some a € [1/2, 1). Note that by letting a = 1/2, 
one retrieves the CIR process. Higham and Mao (2005) study strong convergence in the case of 
(b4). Lord et al.| ( |2010 ) introduce (b2), a modification of (bl), discuss its strong convergence, 
and present an overview of several discretization schemes, including (bl)-(b4) and implicit ones, 
and present numerical comparisons. Similarly, Alfonsi ( |2005 ) presents implicit schemes (that 
admit analytical solutions), studies their weak and strong convergence, and presents numerical 
comparisons with (bl) and (b3). Among popular implicit methods, let us mention the implicit 
Milstein scheme, described for instance in Kahl et al. (2008), where it is found to be better for 
discretizing the CIR process than the explicit Milstein scheme or the balanced implicit method 
of [Milstein et al.| ( |1998[ ). 

In financial engineering, one is often interested in pricing a derivative security for which 
the CIR process is involved in modeling the underlying asset. However, it seems like very 
little attention has been given to the question of convergence of approximate prices (resulting 
from the discretization) to the right price, especially for path-dependent derivatives. In the 
present paper, we address this problem in the even more general framework where the price 
of the underlying asset follows a (typically nonnegative) diffusion process with time- dependent 
coefficients dX{t) = b{t,X{t)) dt + a{t,X{t)) dW{t). The discounted payoff function of, for 
instance, an option is a function of the path of the underlying asset price, say g{X). The price 
of this option is therefore given by E [(7(X)], when the underlying probability measure is a risk- 
neutral measure. The goal is now to define a sequence of approximating processes (Xn)n>i, 
based on a discretization scheme, such that E [(jr(X„)] converges to E [^(X)] as n goes to infinity. 
Indeed, Monte Carlo estimators of E [(/(X„)] will then provide numerical values of the right 
price E[5r(X)]. Thus, existence of a weak solution X to the above SDE, and convergence in 
distribution of a sequence of processes {Xn)n>i to this solution is usually more than sufficient 
for pricing purposes. Indeed, such convergence is equivalent to E [g(X„)] converging to E [^(X)] 
for all path functionals g within a class of sufficiently well behaved functionals. The discounted 
payoff function g of an option is typically a continuous (or almost surely continuous) function 
of the path of the underlying asset, and this is usually sufficient to use the previous definition, 
establishing at once price convergence. It is to be stressed that convergence in distribution of 
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the sequence of processes {Xn)n>i to X involves the distribution of the whole path and must not 
be confused with convergence in distribution of the sequence of random variables (X„(T))„>i 
to X{T), where T is the maturity time. The latter is a much weaker statement and is useful 
for pricing European contingent claims for which the payoff depends only on the value of the 
underlying security at the maturity date, but generally does not allow us to deal with path- 
dependent contingent claims. Weak convergence results generally found in the literature (such 
as in Bossy and Diop ( |2007 ) or Alfonsi (2005)) are of this latter type, i.e., they pertain to the 
processes sampled at a fixed instant T > 0. Note that the convergence in distribution of the 
sequence of processes (X„)„>i to X also includes the convergence in distribution of any sequence 
of random vectors . . . to {X(ti), . . . ,X{tk)) for fixed times ti, . . . ,tk- 

As pointed out earlier, much attention is dedicated to strong convergence in the literature. 
Strong convergence roughly says that the approximating process is uniformly close to X on the 
interval of time [0, T] for large n, hence it holds promises to establish price convergence for path- 
dependent contingent claims with maturity T. Higham and Mao (2005) actually take this next 
step when X is the CIR process; they define continuous-time approximation processes from the 
discrete scheme (b4), prove their strong convergence towards X, and then deduce convergence 
of the price of a few path-dependent derivative securities. Note that many papers deal with the 



CIR process in isolation; see, e.g., Alfonsi (2005), Berkaoui et al. (2008), and Bossy and Diop 



(2007). As opposed to that, Higham and Mao (2005) prove convergence of the price of a barrier 
option in Heston's model, in which the CIR process is used to model the squared volatility 
of the stock price. To the best of our knowledge, they are the first to establish, by showing 
convergence for certain option prices, that using an Euler-type discretization in the full Heston 
model is theoretically correct. Numerical results, obtained from several discretization schemes, 
are provided in Lord et al. (2010) for some options in the Heston model. Note that the transition 
density function of the CIR process is known to be (within a scaling parameter) noncentral chi- 
square, which allows for direct simulation of this process; [Broadie and Kaya ( 2006^ have provided 
an exact simulation algorithm for Heston's model. However, algorithms using the transition 
density are computationally slower than Euler-type schemes, especially when the trajectory 
must be sampled at a large number of time points. Therefore, this family of algorithms is less 
suited for pricing highly path-dependent options; see for example the introductory discussion 



m 



Higham and Mao (2005). For this reason, and the fact that direct discretization methods are 



widely used in practice, our focus is on the latter. 
So 



Higham and Mao (2005) show that strong convergence of the approximating processes 



to X may effectively be used to deduce price convergence for certain contingent claims. Such a 
deduction involves somewhat delicate calculus of probability, the complexity of which depends 
on the complexity of the specific contingent claim considered. As previously discussed, one 
could instead easily deduce price convergence by essentially verifying that the payoff function 
is continuous, provided that one had a sequence of approximation processes (X„)„>i weakly 
converging towards X, in the sense of convergence in distribution of the whole path. It is possible 
to implement this other approach, with relative ease and in general setups, by using a powerful 
idea set forth by Stroock and Varadhan (discussed in detail in Stroock and Varadhan (1979)), 
namely the characterization of Markov processes by means of the so-called martingale problem. 
In particular, Stroock and Varadhan's approach provides an alternative way to regard diffusions. 
This point of view has the advantage of being particularly well suited to establish convergence 
of Markov chains to diffusion processes. Actually, a martingale problem is entirely defined by 
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the expression of a generator, and it turns out that it is sufficient to estabUsh convergence of the 
generators of a sequence of Markov chains to the infinitesimal generator of a diffusion, for this 
sequence of Markov chains to converge in distribution to the diffusion. 

In this paper, our main goal is to propose a discretization scheme, in the form of a Markov 
chain with nonnegative values, and use the above-mentioned techniques based on the generator 
to show its convergence in distribution to the solution of the SDE, under suitable assumptions 
on h and a. Note that in a standard Euler-Maruyama scheme, such as the one presented in ^ 
for the CIR process, it is the use of the normal distribution which is responsible for the non- 
zero probability of getting a negative value on the next time-step, even if the current value 
is nonnegative. When we work within the framework of the martingale problem, there is no 
additional difficulty in establishing convergence if we trade the normal distribution for another 
distribution whose second moment is finite, in the spirit of the weak Euler scheme (see, e.g.. 



Kloeden and Platen (1992)). This suggests the following idea: let us use a scheme very similar 
to the standard Euler-Maruyama scheme, where by a careful choice of distribution we make 
sure that the resulting Markov chain is well defined and assumes only nonnegative values. We 
propose to use a nonnegative distribution, and we give conditions on its mean and variance to 
ensure that the scheme is well defined and converges. Thus, application of our scheme in practice 
reduces to the sole choice of a nonnegative distribution whose mean and variance satisfy a given 
condition, making it relatively simple and versatile. The family of diffusions for which such a 
choice of distribution is possible encompasses several examples of practical interest. For instance, 
our scheme applies to the CIR process, including some cases where the reversion parameter is 



a (possibly correlated) stochastic process (as in Deelstra and Delbaen (1998)), and it can be 
used in the framework of Heston's model. Moreover, verifying that the approximate prices of a 
path-dependent option converge to the real price is then a matter of verifying that the payoff 
functional is continuous (at least on a set of probability 1), and hence one does not have to resort 
to delicate calculus of probability; note that it is left for future work to analytically evaluate the 
rate of weak convergence of this scheme. Finally, it is interesting to note that the convergence of 
the binomial tree towards the geometric Brownian motion, and that of the GARCH(1,1) discrete 
process to the continuous version of this process, are special cases of our scheme. 

The paper is organised as follows. In Section [2| we introduce the nonnegativity preserving 
discretization scheme, and specify for which diffusions it applies. In Section |3j we present the 
main result (Theorem [s]), which says that the approximate processes defined by the scheme 
converge in distribution to the right diffusion process. A brief introduction to the martingale 
problem is provided. In Section |4j we discuss the consequences of the convergence in distribution 
established in the main result. In particular, we see the ease with which it can be used in the 
pricing of derivative securities. Some examples of dynamics for which the scheme applies (and 
converges) are given in Section |5] Finally, the results of two numerical experiments involving 
the proposed scheme as well as schemes (bl)-(b4) are presented in Section |6] 

2. The scheme 
Consider once again the stochastic differential equation (SDE) 
(3) dX{t) = b{t, X{t)) dt + a{t, X{t)) dW{t), t > 0, X{0) = xq, 



where the following is satisfied: 
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Condition 1. The coefficients 6 : M+ x M'^ ^ R"' and a : R+ x ^ M"' are continuous 
functions. Moreover, xq is a constant in M'^. Here, d is a positive integer, M+ := [0,oo), and 
(g) M'^ is the space of d x d-matrices. 

The process (M^(i))t>o represents a d-dimensional standard Brownian motion. The main 
reason to allow equation ^ to be multidimensional is to set a framework general enough to 
include as special cases, for instance, two-factor interest rate models or stochastic volatility 
models, where at least two SDEs are involved simultaneously. We are primarily interested in 
cases where the solution to ([3| is a M'^-valued process, at least one of whose components is a 
nonnegative process. More precisely, we postulate the following: 

Condition 2. There is an integer m, with < m < d, such that for each xq £ E, where 
E := X M'^"™', the SDE ^ has a unique (in the sense of probability law) weak solution, that 
is there exists a probability space (Q, J^, P), a filtration {J-t}t>Q satisfying the usual conditions, 
and a pair [W,X), where W is a d-dimensional standard Brownian motion and X is a process 
with continuous paths satisfying Moreover, the unique weak solution is such that X{t) € E 
for all t > 0, almost surely. 

Remark 3. Note that only the values of the functions b and a over M+ x E are relevant under 
this setting. 

In order to fix ideas, let us first focus on the case where xq > is given, and {X(t))t>o is a 
nonnegative one-dimensional process (i.e., d = m = 1 and E = M+). Moreover, assume that a{-) 
is a nonnegative (scalar) function. If we fix a discrete time-step of size 1/n, for some integer n, 
then a discrete-time approximating process (l^(/s))fc>o is defined as follows: let 1^(0) = xq and, 
for each A; > 0, set 

(4) Yn{k + 1) = Yn{k) + h{k/n,Yn{k)) + -^a{k/n,Yn{k)){e^^ - fi), 

n y n 

where (e[."'*)A:>o is a family of independent copies of a random variable e with mean /i and 
variance 1. In a standard Euler scheme, e — /i follows the standard normal distribution, causing 
Yn{k + 1), even given Yn{k) > 0, to have a non-zero probability of being negative. To avoid this 
problem, let us instead assume that e is a nonnegative random variable. If its (positive) mean 
is set such that 

(5) x+-bit,x) ^cr(t,x)^>0, for all (t, j;) G R+ X R+, 

n \ n 



then (1^(A;))a;>o is clearly a Markov chain whose values are nonnegative. We will see later that 
letting go of the normality is not too much of a price to pay, as (l^(A;))fc>o provides a valid 
approximation (see Theorem [s]). We also delay the illustration of the above simple scheme, in 
the case of the CIR process (in Example [7| , until after its generalization to the multidimensional 
case, that we now undertake. 

Assume the general setting of Conditions [l] and [2] Put a := ao^ (T indicates the transpose 
operation). The multidimensional discretization scheme relies on a function a, a sequence of 
independent copies of a random vector e on a probability space (0, J^, P), and an integer no, 
chosen to satisfy: 
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Condition 4. We have a{t, x) = a{t, x)S {t, x) for all [t, x) E M+ xE, where T, is a symmetric 
semi- definite positive d x d matrix, and a : M+ x E ^ 'K'^ ^ is a continuous function. The 
positive integer uq and the mean n are such that 

(I 1 \ 

(6) inf X H — h (t, x) ]=a it, x) ii] £ E, for all n > no, 

{t,x)m+xE V n yjn ) 



where the infimum is taken componentwise. Moreover, the random vector e has mean fi, covari- 
ance matrix T,, and 

(7) ]P(a(t, x)e G E for all (t, x)eR+x E) = l. 

Remark 5. Equation ^ is a direct extension of condition (jsj). The condition in equation ([7| is 
satisfied, for instance, if all components of a are nonnegative functions, and all components of e 
are nonnegative random variables. When d = 1 and E = and in the typical situation where 
the scalar function a is nonnegative, we usually use a = a and S = 1. In this case, as mentioned 
previously, upon setting a nonnegative distribution for e with variance S = 1 and mean /j, > 
satisfying ^ (or equivalently ^ ), then Condition |^ is satisfied. 

Remark 6. Note that fi and no are typically chosen as a couple in order to satisfy the condition 
given in equation providing more flexibility in the choice of the distribution of e. 

Let (e[,'^'')fc>o, n>no t>e a family of independent copies of the random vector e of Condition [4] 
on the probability space (Q,J^,¥). Fix some xq G E. For n > uq, define (l^(A;))fc>o by letting 
1^(0) = Xo, and then by iterating as follows: 

1 1 

(8) Ynik + 1) = Yn{k) + -b {k/n, Yn{k)) + —a {k/n, Yn{k)) {ef'^ - fi), ke {0, 1,2,...}. 



n ^/n 

For each n > no, the process (5:^n(A;))fc>o is a Markov chain with values in E, i.e., ¥(Yn{k) G 
E) = 1 for all A; > 0, in view of Condition [4] and the fact that xq G E. Finally, define the time- 
continuous approximating process {Xn{t))t>o by linear interpolation of the values of {Yn{k))k>o 
between the discrete time-steps: 

(9) Xnit)=Yn{[nt\) + {nt-[nt\)iYni[nt\+l)-Yni[nt\)), t>0, 

where [y\ is the largest integer less than or equal to y. By construction, {Xn{t))t>o is an i?- valued 
process for each n > no- Our main result, discussed in Section [3j states that the sequence of 
approximating processes (X„)„>„g converges to the (unique) solution of the SDE Section |4] 
discusses further the nature of this convergence and its applications in option pricing. 

We conclude this section by giving an example of dynamics for which Conditions [T| [2] and |4] 
are satisfied, and hence the scheme is applicable. Section [5] gives more such examples of practical 
interest. 

Example 7. Cox-Ingersoll-Ross (CIR) model The CIR process is defined by equation ^ 
with b{t,x) = k(/3 — x) and a{t,x) = V\fx, x > 0, where k, /3 and v are positive constants. 
Here, d = m = 1 and E = M_|_, and we trivially set a = a and S = 1. With straightforward 
optimization computations, one can show that whenever n > k and n > we have 

(10) inf fx + -K(/3-x)-^z.VS^ ^'"^ 



x>Q \ n -y/n J n 4(n 
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Consequently, if we choose no > k, and /x such that the right-hand side of ( 10 ) is nonnegative, 

and finally set a nonnegative distribution with mean fi 



which is true ifO<fi< k(3 (^1 — - 

and variance 1 for e, then Condition^is satisfied. It is also interesting to note that Js]) implies 
E[Yn{k + 1)] = nP/n + E[Yn{k)]{l - n/n), = 0,1,2,..., from which we deduce E[y„(A;)] = 
/? + (xo - /3)(1 - K/nf. Using one may easily conclude that, for a fixed time T > 0, 
lim„^oo]E[X„(r)] = (3 + {XQ - f3)e- ''^ = E[XiT)], where {J^t}t>o, {W,X) is the weak 

solution of ([3]); see (Shrev^ 2004b Section 4-4) for the last equality. Very similarly, one may 
establish that the variance of Xn{T) goes to that of X(T) as n goes to infinity. 

3. Convergence of the scheme using the martingale problem formulation 

In this section, we establish convergence of the probabihty law of the processes X„, n > no 
(defined in the previous section), towards the law of the solution to the SDE ([3|, as n goes to 
infinity (or the time-step goes to zero). We achieve this by means of the martingale problem of 
Stroock and Varadhan. We provide a very brief introduction to the martingale problem in this 
section, mainly establishing the notation, and refer any reader seeking for a detailed discussion 
to Ethier and Kurtz (1986) or Stroock and Varadhan (1979). 

Let C^{R'^) be the set of infinitely differentiable functions / : R"' — > IE 
Define the differential operator A, acting on functions / G C^(M'^), by 



with compact support. 



d ^ d d 

(11) {Af){t,x) := j;6,(i,x)5,J(x) + 2 j;^a,j(t,x)9,,a,^,/(x), 

i=l i=l j=l 

(t, x) G M+ X M'^, where dx^ stands for the partial derivative with respect to the i-th variable, 
and hi and aij denote the entries of b and a. If (Jl, P), {J't}t>o, iW,X) is a weak solution to 
the SDE ([3]), then P(X(0) = xq) = 1, and from Ito's formula one easily deduces that 

(12) f{X{t))- [' Af{s,X{s))ds 



is an {J^f}-martingale for any / G C^(M'^). Actually, any process X with continuous paths, 
defined on a probability space (fi, J^, P) endowed with a filtration {J-^t}t>o, satisfying P(X(0) = 
Xo) = 1 and such that (12) is an {J^^j-martingale for all / G C^(M'^), is called a solution 
to the martingale problem for (A,xo). Hence, any weak solution to the SDE ([3| solves the 
martingale problem for {A,xo). Under Condition [T] the converse turns out to be true. In fact, 
if the martingale problem for {A, xq) has a solution, then there exists a weak solution to the 
SDE ^ (see Corollary 5.3.4 in [Ethier and Kurtz| fl986| )). Moreover, uniqueness (in the sense 
of probability law) holds for solutions of the SDE (|3|) if and only if uniqueness holds for the 
martingale problem for (^,xo) (see again Corollary 5.3.4 in Ethier and Kurtz (1986)). In other 
words, existence and uniqueness of a weak solution to the SDE (|3|) is equivalent to existence 
and uniqueness of a solution to the martingale problem for {A,xo). This means that all the 
information which is truly essential relatively to the SDE ^ is encapsulated in the martingale 
problem formulation, which in turn relies only on the expression of A, also called the infinitesimal 
generator. 

In view of the importance of the martingale problem and the generator, let us expose similar 
ideas in discrete time, that is, let us see how the Markov chain (l^(A;))fc>o defined in (^ may be 
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characterized through a martingale problem as well. Consider the transition function 

(13) KJt, X, r) := IP ( X + ^b(t, x) + -^a(t, x)(e - fi) eV] , T G BCM.'^), 

V n \ n I 



for (t, x) G M+ X where i3(M'^) represents the Borel sets on W^. This function actually describes 



the transitions of the Markov chain (1^(/c))a;>o- Indeed, let {-7-"^"}fc>o be the filtration generated 



by (l"n(fc))fc>o (i-e. := cT{y„(/)|/ = 0, 1, . . . , /c}), and note that 

(14) ]P(y„(A; + 1) G r I Tl-) = Knik/n, y„(A:), T), T G B{R'^), 

for G {0, 1, . . . }. For all / G C^°°(M"'), define 



Anfit, x):=n {f{y) - f{x))Kn{t, x, dy), 



(15) 



(t, x) G M+ X E. For each / G C;?^(M'^), it is easily seen that 



(16) 



fc-i 

fiYn{k))--Y,Anf{l/n,Ynil)), 



1=0 



is a {J^^"}-martingale. Conversely, any discrete-time process (Yn{k))k>o such that Yn{k) = xq 
and ( 16 ) (with Yn in place of 1^) is a martingale for all / G (MJ^) (with respect to the filtration 



generated by {Yn{k))k>o) is a Markov chain whose transitions are described by Kn (Stroock 



and Varadhan 11979 Section 11.2). Hence, the martingale problem defined by (15) and (16) 



characterizes such Markov chains, and An in (15) is the discrete analogue of the infinitesimal 
generator A. 

In order to show convergence of the sequence of processes {Xn)n>no (recall (|8| and (|9l)) to 



the weak solution of the SDE (|3]) as n goes to infinity, it turns out that it is sufficient to show 
convergence of the sequence of generators, more precisely that A^f converges to Af as n goes 
to infinity, for all / G C^iW^). It is this method that we use to establish our main result: 

Theorem 8. Under Conditions^ [£| anc?[^ the sequence of approximating processes {Xn)n>no 
defined by ([8| and ([9| converges in distribution to the weak solution of SDE ^ or, equivalently, 
to the solution of the martingale problem for (^, xq)- 

The details of the proof are found in the appendix. We insist on the fact that the convergence 
in distribution established in Theorem [8] is that of the whole path of the processes X„, n > no, 
to the whole path of the solution to the SDE ([3|, although we delay until the next section a 
careful definition of this convergence. 



4. Consequences of the convergence in distribution 

Let us take the example where the SDE ^ models the evolution of stock prices in the risk- 
neutral world. Let (J7,J^, P), {J^t}t>Oi {W,X) be the weak solution of ([3|. Suppose that the 
discounted payoff of a derivative security is expressed as a function, say h{X), of the underlying 
price process X] note that path-dependent options are embedded in this setup. On one hand, 
it is well known that E[/i(X)] is the price of this derivative. On the other hand, saying, as 
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in Theorem [sj that the sequence of approximating processes {Xn)n>no (defined on {(},J^,F)) 
converges in distribution to X means by definition that 

(17) limE[g{Xn)]=E[g{X)], 

n— ^-oo 

for all bounded continuous functions g. Thus, if h is nice enough (this is true in particular if h 



is itself bounded and continuous), then (17) holds with g = h, and the approximate price given 
by the scheme, i.e., E[/i(X„)], is close to the actual theoretical price when n is large. 

For the definition of the convergence of {Xn)n>no to X to be complete and transparent, this 
section starts by clarifying what it means for the above-mentioned functionals g to be continuous. 
Then, it provides sufficient conditions for h to be nice enough for the sequence of approximate 
prices to converge to the right price. Moreover, this section includes several examples of such 
nice functions, and illustrates how they can be used for pricing. 

The solution of the SDE ([3| and approximating processes X^ constructed in Section |2] are all 
processes whose paths are continuous functions. In other words, these paths belong to CeO^+), 
the set of continuous functions x : M+ — )• E. Let us endow Ce(^+) with the topology of uniform 
convergence on compact subsets of M_(_, induced by the metric 



Aix,y) := / e-" sup {\x{t) - y{t)\ A 1) du , 
Jo 0<t<u 

for X, y E C£;(M+) (where | • | is the Euclidean norm on M'^ and uAv is the minimum of n, f E M). 
The processes Xn, n > no, and X may be regarded as random variables with values in the set 



C£;(M+), and convergence in distribution of X^ to X means that (17) is satisfied for all functions 



g : Ce(^+) — K which are bounded, and continuous with respect to the metric A. 

Remark 9. Consider a positive constant T and a function g : CeO^+) — >• IS. // g{x) depends 
only on the values x{t), <t <T, for any x E CeO^+), then g is continuous with respect to the 
metric A if and only if it is continuous with respect to the usual and more tractable sup-metric 
defined as 

Ar(x,y):= sup \x{t)-y{t)\, 

0<t<T 

for x,y E Ce{^+)- In mathematical finance, when we consider a finite investment horizon [0,T], 
the discounted payoff of a derivative security is usually of the form g{X), for some function 
g : Ce{^+) — M whose values depend only on the trajectories restricted to the interval [0,T], 
and hence one simply has to verify continuity of g with respect to At- 



We can enlarge the set of functions g such that (17) holds, extending at once the set of admis- 



sible payoff functions. Indeed, as a direct consequence of Theorem ^ together with Theorems 



1.5.1 and 1.5.4 in Billingsley (1968), we get: 



Proposition 10. Assume Conditions [i| [£| and Let (X„)„>„(, be the sequence of processes 
defined by ([s]) and ([9| on the probability space {Q,T,'F), and let {Q,T,F), {J-t}t>0) {W,X) be 
the weak solution to the SDE ([3|. Suppose that g : Ce{^+) — )• M is continuous with respect to A 
except on a subset Dg C CeO^+) satisfying ¥{X E "Dg) = 0. // the sequence of random variables 
{g{Xn))n>no is Uniformly integrable, then '&\g{Xn)] — ?• E[5f(X)] as n ^ oo. In particular, this is 
true when g is bounded. 
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Remark 11. For a given n, the central limit theorem (or the strong law of large numbers) 
implies that the average of the payoff g{Xn) of N Monte Carlo simulations of the trajectories of 
Xn converges to K[g{Xn)]- The latter is in turn very close to K[g(X)] when n is large in view 



of Proposition 10 So, the price of a derivative security is obtained in practice by letting both 
n and N tend to infinity. This is very similar to the binomial tree method used in the Black- 
Scholes model ( whose connection with the proposed method will actually be discussed further in 



Section 5.1) 



For a positive constant T (typically the investment time horizon), a few examples of contin- 
uous payoff functionals g : C£;(M+) — )• M are given by g{x) = xi(T), g{x) = (^j^ xi{t)dt^ /T, 

and g{x) = maxo<f<T where xi is the first component of x. Their continuity with 

respect to At is easily verified. These functions are useful when dealing respectively with 
plain vanilla European, Asian, and lookback options. Let us further consider the functions 
g{x) = I{maxo<t<T a^i(i) < -B}, where B is a constant and I denotes the indicator function, 
which is continuous except on the set T>g = {x S Ce{^+) ■ maxo<t<T a^i(i) = B}, and 
g{x) = I{xi(T) < B}, continuous except on T>g = {x E Ce(^+) ■ xi{T) = B}. These are 
typical of situations where barrier and binary options are involved. Finally, note that the map- 
ping g : Ce(J^+) — C]r(M+) defined by g{x) = e^^ for x £ C^(M+) is continuous. This last 
mapping may be useful when the log-price of an asset is modeled by the SDE ([3|. Combining 
functions such as the above allows one to express a wide variety of bounded payoffs, and then 



easily conclude by Proposition 10 that the sequence of approximate prices for the correspond- 
ing derivative security converges to the right price as the time-step goes to zero. Example [12] 
illustrates this technique. 

Example 12. Consider the function g : Ce(^+) — )• K defined by 

9{x) = exp XI (s) ds^ x(^K-^£ e"^^^) ds^ , 

for X £ Ce(^+), where K and T are positive constants and Xi is the i-th component of x. It is 
continuous over Ce(J^+), and bounded if m > 1 (i.e. xi{t) > for all t > 0). Thus, it satisfies 



the assumptions of Proposition 10 If the first two components of a multi- dimensional process 
X, namely Xi and X2, correspond respectively to the short-rate and the log-price of an asset, 
then g{X) is the discounted payoff of an Asian put option with strike price K on this asset. 



We can similarly tackle all the examples considered in Higham and Mao (2005) since these 
are for path-dependent options with bounded payoffs, for instance up-and-out call options. Un- 
bounded payoffs can be treated via the put-call parity principle or by truncation, as in, respec- 



tively. Examples 03^ and 14 below. 



Example 13. Assume Conditions'^ and^ with m > 2. Suppose that ^ models the risk- 
neutral evolution of X, whose first component Xi is the (nonnegative) short-rate and second 
component X2 is an asset price. Let D{t) = exp(— X\[s)ds), t > 0, be the discount factor 
process. In the postulated risk-neutral world, {D{t)X2{t)}t>o is a martingale. Upon combining 
this with equality X2(T) - K = {X2{T) - K)+ - {K - X2(T)) + , we get 

(18) E[Z)(r) {X2{T) - K)+] = X2(0) - K¥.[D{T)] + E[D{T) {K - X2{T))+l 



12 



A SIMPLE DISCRETIZATION SCHEME 



where T > is the investment horizon and K > the strike price. On the right-hand side, both 



payoffs are bounded and continuous, hence we can apply Proposition 10 to approximate the prices 
of the bond and the put option. The price of the call option is then deduced from the put-call 
parity (18). Furthermore, if ¥{maxQ<t<T ^2{t) = B) = 0, where B > 0, then, for example, one 



can similarly price an up-and-in call barrier option from the call option price evaluated in (18) 
and the price of an up-and-out call option whose payoff is bounded: 



E 



D{T) {X2{T)-K)+l 



max X2{t) > B 

0<t<T 



E[D{T) {X2{T)-K)- 



■E 



D{T) {X2{T) - K)- 



li max X2{t) < B 
|^0<t<T 



Example 14. For simplicity, assume that d = 1. For K > 0, let h{z) = [z — K)^ , z G M, and 
let (/ifc)fc>i be an increasing sequence of nonnegative, continuous and bounded functions on M 
converging pointwise to h; for example, one can take hk{z) = [z — K)^ A k, k = 1,2, . . . . As 
TTTix) = x{T) is a continuous function over Ce(^+), then g = ho ttt and Qk = h^o ttt oltg also 
continuous functions over Ce{^+), o-nd in the latter case it is also a bounded function. Then, 
by Proposition\l 0\ for each k we have 

lim E[hk{Xn{T))] =E[hk{X{T))]. 

n— >oo 

Finally, by monotone convergence, E[hk{X{T))] tends to E[h{X (T))] as k goes to infinity. Con- 
sequently, for k and n large, E[hk{Xn{T))\ is arbitrarily close to the price of the call option with 
maturity T and strike price K. 

Remark 15. In order to establish convergence of the sequence of approximate prices to the 
right price for a derivative security (with a bounded payoff function), an alternative strategy is 
to first show strong convergence of the sequence of (continuous-time) approximate processes to 
the solution of the SDE then find a way to deduce, from this strong convergence result, 
price convergence. This is the approach used by Higham and Mao (2005) when ^ corresponds 
to the CIR dynamics. As seen in Higham and Mao (2005), the second step (the one in which 
price convergence is deduced) involves probabilistic arguments, which are specific to the given 
derivative security, and whose complexity vary depending on that of the given derivative security. 
One source of difficulty comes from the way the discrete-time approximation is completed in 
between the discrete time steps to get a continuous-time process strongly converging to the CIR 
process. Indeed, this completion relies on a Brownian trajectory (see (16) in Higham and Mao 



[2005), and strong convergence results Theorems 3.1 and 3.2 and Corollaries 3.1 and 3.2). Such 
trajectory is known only in theory, while practical computation of the Monte Carlo approximate 
price requires a trajectory known in practice. As a result, the trajectory used for the strong 
convergence results is different from that used in the expression of derivative 's Monte Carlo price 
(compare (16) with (18) and, for instance, (26) or (28) in Higham and Mao ^2005 )). Hence, 
the analysis necessarily requires establishing that these two trajectories are close (see Lemma 3.2 
in \Higham and Mao (2005)). In contrast, in the weak convergence approach proposed here, 
the continuous-time approximate processes whose convergence to the solution of ^ is shown 
is known in practice: it is obtained by simple linear interpolation from the discrete scheme 
(recall ^). More importantly, upon availability of such a weakly convergent process, one gets 
price convergence essentially by verifying that the payoff function is a continuous function of the 
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path. As we have shown, this is quite simple and technical probabilistic arguments typical of the 
strong convergence approach are conveniently avoided. 



5. Examples of models and diffusions 

We now give several examples of models and diffusions for which Theorem |8] applies and hence 
the scheme converges. 

5.1. Black-Scholes model with time-dependent coefficients. Let /3(t) and z/(t), t > 0, be 
continuous functions, with u(t) > for all t > 0, and suppose that xq > is given. Assume that 
sup^>QZ/(t) G (0,oo), and mit>o(3{t) > -co. Upon letting b{t,x) = x/3{t) and a{t,x) = xu{t), 
then the solution to the stochastic differential equation ^ is the one-dimensional geometric 
Brownian motion (here, d = m = 1). Set a = a, and S = 1. Then Condition |4] is satisfied if e is 
a nonnegative random variable, and no and /i > are chosen such that 

x(l + — - > for alH > 0, X > 0, n > no . 

V n \/n I 



This will be the case if, for instance, no > — 2 inft>o /3(t) and < < y^/(2 sup^o i^(t)). For 
each n > no, equation ^ gives the following scheme to approximate the geometric Brownian 
motion with time-dependent coefficients and with initial value xq: set 5^(0) = xq and then, for 
A; > 0, set 

(19) Y^ik + 1) = Y^ik) f 1 + ^ + ^^(4") - ,) 



n A/n 



One possible choice of distribution for e is 
(20) F{e = 0) = ^ and F(e = fi+ * 



or, in other words, fie/{fi^ -|- 1) ~ Bernoulli(/i^/(l -|- ^^)). One can verify that we have indeed 
E[e] = ^ and Var(e) = 1 (= S). If the coefficients are constant, i.e. f3{t) = Pq G M, and 



iy{t) = uq > 0, for all t > 0, and e follows the distribution in (20), then it is interesting to note 
that the scheme ( |19[ ) reduces to a recombining binomial tree, in which Yn{k + 1) is equal to either 
UnYn{k) or dnYn{k), where u„ = 1 ^ and dn = I + ^ - ^/U. By construction, we 

have < dn < 1 + ^ < Un, and it is easy to verify that 

,,2 (\ A_ §0) _ r] 

(21) FiYn{k + l) = uMk)) = = ^ 

1 + Un-dn 

If Pq stands for the constant interest rate, and {l^(A;)}fc>o models the risk-neutral evolution of 
an asset price, then the probability of an up-move (equivalent to multiplying the current price 



by Un) found in (21 ) is consistent with the risk-neutral probability usually postulated in binomial 
models (see for instance Shreve ( 2004a| )). 



14 



A SIMPLE DISCRETIZATION SCHEME 



5.2. Constant elasticity of variance (CEV) model. It is possible to apply Theorem [s] to 
one dimensional diffusion processes as those in Berkaoui et al. (2008) which have a non-Lipschitz 



diffusion coefficient. Let b{t, x) = b{x) be a Lipschitz continuous function such that 6(0) > and 
let a{t,x) = ux", with i/ > and a € [1/2, 1). Under these assumptions, and with xq > 0, there 
exists a nonnegative strong solution to equation ^ ( [Berkaoui et al. 2008), and so d = m = 1 
and E = M+. Let e be a nonnegative random variable with mean // and variance S = 1, and 
take a = a. It remains to show that no and fi > can be chosen such that the condition in 
equation ([6| is satisfied, namely 



(22) 



b(x) , , 
— + c„ X > 0, 
n 



for all a; > and n > no, where c„(x) := x — 

\b{x) — b{y)\ < K\x — y\ for all x,y £ M), then b{x) > 
immediate if —Kx/n + Cn{x) > 0, which is equivalent to 



fi. If b has Lipschitz constant K (i.e. 



-Kx, for all x > 0. Thus, (22) is 



n 



(1 - K/n) 



-1 



1/(1-") 



when 1 — K/n > 1 — K/no > 0. Now, it remains to secure inequality (22) for x G [0, x^"-*). 

Since x^"^ — )• when n — )• oo, and since b is continuous, then for any ^ we can set no large 
enough to ensure both no > K and 



(23) 



min b{x) > ^'^^^ n > no- 

3;e[0,a:(")] 2 



We have c'„(x) = 1 - ^"i-c fJ- and c"(x) 
with a global minimum in 



'^'r^\-} IJ- Therefore, the Cn's are convex functions 



/ \ l/(l-a) 



with value 
(24) 



\ l/(l-a) 

n ) 



(^aV(i-)_ <o. 



where the last inequality follows from a < 1. When a G (1/2,1), we have |nc„(x„)| — )• as 
n — )• cxD, and for any ;U one may choose no sufficiently large to get |nc„(x„)| < 6(0)/2 for all 



n > no. Together with (23), this gives 
b(x 



(25) 



+ c„(x)>- mm 6(x) + c„(x„) > - — - +nc„(x„j > 0. 
n ' ' n xe[0,x(")] ^ ' n \ 2 ' ' ' 



for all X G [0, x^")] and n > no. When a = 1/2, then c„(x„) = -i/2/i2/(4n) (recall ([24|). Then 
the last inequality in (25) is immediate if fx satisfies < y'2b{Q)/v. 
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5.3. One-dimensional affine diffusion process. Let ho,hi,ko,ki,rQ be constant numbers 
satisfying kohi — kiho > 0, hi ^ 0, and ho + hiro > 0. Then, the process defined by 

dR{t) = {ko + kiR{t)) dt + ^ho + hiR{t)dW{t), R{0) = ro, 

belongs to the family of affine diffusion processes (see Duffie et al.| (2000)) and takes its values in 
[— /io//ii, oo) if /ii > 0, or in (— oo, — /io//ii] if /ii < 0. The process (X(t))i>o defined by the change 
of variable X{t) = ho+hiR{t) satisfies the dynamics ^ with b{t, x) = b{x) = {kohi — kiho)-\-kix, 

with a = 1/2, 



5.2 



a{t,x) = \hi\^/x, and xo = ho + hirQ. This is the special case of Example 
6(0) = kohl — kiho > 0, i' = \hi\, and Lipschitz constant K = \ki\. In view of Example 5.2 
we must choose /x < y^26(0) /u = y^2{kohi — kiho)/\hi\, and then take no > K large enough 
for (23) to be satisfied, for instance no > max(2Er, 8|/ci|i^^/U^/6(0)). 



5.4. Two-factor CIR model. In a two-factor interest rate model, the interest rate process is 
defined as r{t) = 5o + 8iXi{t) + 52X2{t), where (5o > and (5i, (^2 > (see|Shreve| (|2004b[)). In 



b{x) 



the two-factor canonical CIR interest rate model, the two-dimensional process {X{t))t>o evolves 
according to the general dynamics given in equation ([3|, where the coefficients are 

/3i - AiiXi + A12X2 \ ^ / y/xl 

P2 + A21X1 - A22a:;2 J V v^\/l ~ 

for X = (xi, X2) £ E = M^, and with initial condition X{0) in M^. The parameters /3i, P2, An, A22 
are positive constants, A12, A21 are nonnegative, and the instantaneous correlation p satisfies 
— 1 < p < 1. Note that we are in the case d = m = 2. We define, for x G E, 

^/x^ J \ P ^ 

Very similarly as in Example [TJ we may conclude that Condition [4] is satisfied as soon as 



(26) 



a{x) 



(27) 



no > max(Aii, A22) and < m < 2 




1,2, 



and the 2x1 random vector e has nonnegative components. Note that if e is a nonnegative 
random vector with E[e] = /i, and covariance matrix S specified in (26), then it must be true 
that 



(28) 



- o< 



Mim- 



Conversely, under condition (28), one may construct a random vector e with nonnegative com- 
ponents, E[e] = fi, and covariance matrix S (specified in (26)). As a result, it is possible to 
choose no and ^1 , /i2 such that (27) and Condition 4 are satisfied as long as —p < A\/ fiiP2- 



5.5. Stochastic volatility models. Let a,A,i/ be positive, /3 G M, and p G (—1,1). Consider 
the G ARCH (1,1) stochastic volatility model 



(29) 



dV{t) = {a- \V{t))dt + vV{t)dWi{t); 
dS{t) = S(t){pdt + ^/v{i)dW2{t)), 



with ^(0) = vq > Q and 5'(0) = So > 0, and where Wi and W2 are two standard Brownian 
motions with instantaneous correlation p. Even though both processes are nonnegative, we ease 
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the situation by working (as in Lord et al. (2010)) with log 5" instead of S, which does not need 



to remain positive. More precisely, we consider equation ^ with d = 2 but m = 1, and 



for X = (xi,X2)^, and with X{Qi) = (uq, log(so))^- Define 

(30) *,.).(7 » ) and ^-{ll 



Note that ([7| in Condition |4] holds as soon as ei is a nonnegative random variable (e2 may take 
both positive and negative values). For Condition |4] to be satisfied, it remains to make sure that 

inf ( xi H ^ 7=/^i ) ^ ) for n > no • 

xi>a \ n J 

This will be the case if we choose ng > A and < ;Ui < ^1 — . Note that, for n > no, 
the approximating scheme for {Xi{t))t>o = {V{t))t>o given by ^ may be written as 

+ 1) = ^ + y^Ak) (i - ^ - ^/ii) + ^yuk)e^;:l , 



(where Yn = (yi,n,^2,n)^ and e^""^ = (ei"fc\ £2"^)''^) ^^"^ ^i,n(0) = vq. This is commonly refered 
to as a GARCH(1,1) discrete process, with the particularity that the family of random variables 
i^i^k^k>o are i.i.d. and follow any nonnegative distribution with mean /ii and unit variance. 
Hence, convergence of the GARCH(1,1) discrete process to its continuous counterpart is a special 
case of Theorem [H 



Note that if the process in (29) were replaced with a CIR process, a very similar argument 



could be applied, which would lead to a scheme for Heston's stochastic volatility model. 



6. Numerical results 

In this section, we present the results of two numerical experiments putting the proposed 
discretization scheme to the test. Results pertaining to other discretization schemes are also 
provided in order to facilitate comparisons. In the first experiment, we shall be interested in 
the price of a bond in the CIR interest rate model. This is a good test case since there is 
a well-known closed form formula for the bond price allowing for bias computations, and the 
payoff is path-dependent allowing to test the ability of the method for pricing path-dependent 
derivatives. In the second experiment, we consider the price of a (plain vanilla) European call 
option in Heston's stochastic volatility model. This is another interesting test case as it involves 
the joint action of two correlated processes, the CIR process being one of them, while the pricing 



problem is still mathematically tractable enough to obtain theoretical prices (Heston (1993) 
determines the characteristic function of the stock log-price at maturity; Fourier inversion can 
thereafter be used to get the price of a European call option). 
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6.1. Bond pricing in the CIR model. For the first experiment, we consider tlie CIR process of 
Example|7j tliat is {X{t))t>o is a solution of ^ wlien b{t, x) = K,{f3 — x) and a{t, x) = v^/x, wliere 
K, P and 1/ are positive constants. Assuming tliat {X{t))t>o models the interest rate, the value of a 



exp ^— X{s) ds 



bond that will pay its face value of 1000 dollars in T years is worth 1000 xE 

today. Our goal is to get an approximation of this price using Monte Carlo simulations together 
with various discretization schemes and compare them to the true value of the bond; see, e.g., 
( [Lamberton and Lapeyre 1991, Proposition 6.2.5) for the closed-form formula. 

It is well known that the CIR process will never hit the origin, that is ¥'{X{t) > for all t > 
0) = 1 if xq > 0, when the model parameters satisfy \/2kI3 > z^, while it eventually hits 



the origin with probability one when -v/2k/3 < i^; see, e.g., (Lamberton and Lapeyre 1991 



Proposition 6.2.4). Since we want to assess how the different discretization schemes manage the 
zero boundary condition, it is natural to test the methods for parameters for which the boundary 
problem will theoretically arise. In fact, we put k = 0.5 and /3 = 0.04, so that ^2k/3 = 0.2, 
and consider two values exceeding 0.2 for the volatility coefficient: v = 0.3 and v = 1. In the 
latter case, the zero boundary problem is magnified by the larger volatility and so it should be 
intuitively more difficult to get accurate prices. Tables [T] and [2] show, for different values of n 
and Euler-type discretization schemes, the bias and the margin of error at the 95% confidence 
level for the Monte Carlo prices, based on one million trajectories of the CIR process (i.e., 
N = 1 000 000 in the notation of Remark 11 ). For the reader's convenience, those biases that are 
not significantly different from zero at the 95% level appear in bold fonts. In our discretization 
scheme (the one defined in (|4])), we have set for £ the two- valued distribution, one of whose 
values is zero, with average ^ = 0.8 in the low volatility case and ji = 0.28 in the high volatility 
case (recall (20), where this Bernoulli-type distribution is explicitly given). Note that with 
these values of /i the condition < < ^ ('^/3(1 — ('^/^o)))^''^ from Example [t] is satisfied with 
respectively no = 4 and no = 50 for the two sets of parameters. This choice of distribution is 
motivated by the fact that the Bernoulli distribution is the simplest there is, which results in 
reduced simulation time. Note that the approximation of the payoff term X{s)ds requires 
an approximation of the path based on the discretization scheme for the entire interval [0,T]. 
For the scheme (b4), the trajectories have been interpolated in between the time-steps using 



Xn{s) = \Yni[ns\)\, s G [0,T], since (Higham and Mao 2005 (4), (26), Theorem 4.1) prove 



that the bond approximate prices converge to the right price as n goes to infinity under this 
setup. For all other schemes, the linear interpolation defined in ([9| is used. Since the payoff 
function g{x) = 1000 x exp{— x{s) ds), for x £ Cir_|_(M+), is continuous and bounded, then 
Proposition [To] implies that Monte Carlo bond prices from our scheme converge to the right price 
as n goes to infinity. 

In Table [T] and Table [2] we see that biases are consistently less than a dollar for our scheme 
(in the Bernoulli column) and schemes (bl) and (b2), and for the values of n shown; this is 
quite small in comparison to the true bond prices which are respectively 925 and 940 dollars. 
The biases and intervals at the 95% confidence level of those three schemes are then compared 
graphically in Figure [T] and Figure |2] As anticipated, for all methods it is more difficult to 
evaluate the bond price in the high volatility case. Indeed, in the latter case much more time 
steps per year are required to get an approximate price which is not significantly different from 
the true price (at the 95% confidence level, based on one million trajectories). 
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Table 1. Bond pricing in the low volatility case 
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(U.lio ) 


(U.llzj 


10 


U. lO-y: 


n 1 HQ 


U.ao / 


— D.U / U 




(0.120) 


(0.121) 


(0.122) 


(0.115) 


(0.114) 


20 


0.060 


0.057 


0.166 


-3.545 


-1.860 


(0.119) 


(0.120) 


(0.120) 


(0.116) 


(0.116) 


40 


-0.062* 


0.057* 


-0.051 


-2.257 


-0.971 


(0.119) 


(0.119) 


(0.120) 


(0.117) 


(0.117) 


80 


0.029* 


-0.008* 


-0.082* 


-1.558 


-0.373 


(0.118) 


(0.119) 


(0.119) 


(0.117) 


(0.118) 


160 


0.008* 


-0.000* 


0.005* 


-0.991 


-0.261 


(0.119) 


(0.119) 


(0.119) 


(0.118) 


(0.118) 


rate of conv. 


0.54 


0.94 


1.41 


1.99 


2.85 



Bias and margin of error at the 95% confidence level I 
CIR process parameters: k = 0.5, l3 = xq = 0.04, v = 
Bond parameters: T = 2 years and face value of 1000, 
True bond value: 925.258. 
Mean for e: n = 0.8. 



in parent 
0.3. 



dieses). 



Table [T] and Table |2] were designed primarily to compare biases, some of which are quite 
different in sizes for the same n. As a complement of information, their last rows provide a rough 
estimate of the weak order of convergence. Recall that these tables are based on one million 
trajectories, which is actually not sufficient to get precise numerical orders of convergence. Here, 
the rate of convergence is estimated to be the slope (in absolute value) , when linearly regressing 
log(|bias|) on log(n) (and a constant term). For methods that do very well, note that several 
values of the bias shown in the tables are not significantly different from zero. Hence, these 
values carry very little information on the exact size of the deviation between E[5'(X„)] and 
the real price E[5((X)]; they are essentially noise. Hence, some of these values (marked by 
an asterisk in the tables) have not been taken into consideration in the regressions. By using 
much larger numbers of trajectories, the actual deviation between E[(5r(X„)] and the real price 
has been estimated with much more accuracy in the case of the proposed scheme. The results 
are reported in Table [3j In the low volatility case, the estimated order of convergence of 1.2 
lies between the rough estimates obtained for methods (bl) and (b2), namely 0.94 and 1.41. 
However, the estimated order of convergence of 0.36 is smaller than both the rough estimates 
for methods (bl) and (b2) in the high volatility case, so the actual order of convergence might 
be slower for the proposed method. Calculation of the real order of weak convergence could be 
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Figure 1. Bias comparison in the low volatility case 



1.6 
1.4 
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1/n 
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Table 2. Bond pricing in the high volatility case 



n = steps/year 


Bernoulli 


(bl) 


(b2) 


(b3) 


(b4) 


50 


-0.678 
(0.249) 


2.044 
(0.270) 


4.720 
(0.271) 


-117.019 
(0.318) 


-108.046 
(0.311) 


100 


-0.329 
(0.250) 


0.798 
(0.263) 


2.086 
(0.263) 


-98.960 
(0.313) 


-90.511 
(0.307) 


200 


-0.435 
(0.252) 


0.442 
(0.259) 


1.263 
(0.257) 


-84.907 
(0.308) 


-76.459 
(0.302) 


400 


-0.368 
(0.253) 


0.278 
(0.257) 


0.453 
(0.257) 


-74.073 
(0.304) 


-66.008 
(0.298) 


800 


-0.207 

(0.253) 


0.007 

(0.256) 


0.144 

(0.256) 


-65.310 
(0.300) 


-57.324 
(0.294) 


1600 


-0.064* 

(0.253) 


-0.067* 

(0.256) 


0.402 
(0.255) 


-58.399 
(0.297) 


-50.397 
(0.290) 


rate of conv. 


0.33 


1.77 


0.88 


0.20 


0.22 



Bias and margin of error at the 95% confidence level (in parentheses). 
CIR process parameters: k = 0.5, l3 = xq = 0.04, v = 1. 
Bond parameters: T = 2 years and face value of 1000. 
True bond value: 940.024. 
Mean for e: ^ = 0.28. 



the subject of future work. Note that regression of the root mean square error (RMSE) instead 
of the log-bias would have produced similar results. 
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Figure 2. Bias comparison in the high volatihty case 



2 




0.005 



0.01 0.015 
1/n 
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Table 3. Bond pricing using the Bernoulh distribution 



low volatihty case 




high volatility 


case 




n 


N 


bias 


margin 


RMSE 


n 


iV 


bias 


margin 


RMSE 


4 


4 X 10^ 


0.1951 


0.0616 


0.1977 


50 


4 X 10*^ 


-0.4800 


0.1243 


0.4842 


8 


16 X 10^ 


0.0675 


0.0302 


0.0692 


100 


8 X 10*^ 


-0.3940 


0.0884 


0.3965 


10 


25 X 10^ 


0.0801 


0.0240 


0.0810 


200 


16 X 10*^ 


-0.2883 


0.0628 


0.2901 


20 


100 X 10^ 


0.0197 


0.0119 


0.0206 


400 


32 X 10^ 


-0.2583 


0.0447 


0.2593 


40 


400 X 10^ 


0.0134 


0.0059 


0.0137 


800 


64 X 10^ 


-0.1694 


0.0317 


0.1701 


rate of conv. 


1.209 




1.201 




0.362 




0.363 



Bias, margin of error at the 95% confidence level, and RMSE for different combinations 

of n (the number of steps per year) and (the number of trajectories). 

The parameters are the same as in Table [T] (resp. Table |2]) in the low (resp. high) volatility 

case. 



From this first numerical experiment, one can safely conclude that our scheme is very compet- 
itive when it comes to discretizing the CIR dynamic and pricing a path-dependent derivative, in 
both low and high volatility environments. 

6.2. Pricing of a European call in Heston's model. For the second experiment, we consider 
Heston's stochastic volatility model (in the risk-neutral world), in which the squared volatility V 
and stock price S evolve according to 

dV{t) = K{^-V{t))dt + uy^V{ijdWi{t), 

dS{t) = S{t){rdt + y^dW2{t)), 
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with V{0) = vq and S{0) = sq. Here, k, /3, z^, r, and sq are positive constants, and Wi 
and W2 are two standard Brownian motions with instantaneous correlation p. Our goal is to 
approximate the price of a European plain vanilla option whose value is K[e~'^'^ {S {T) — K)~^]. 
For ease of comparison, we use the same set of model parameters as in experiment SV-I in [Lord 
et al. ( 2010[ ) and experiment two in Broadie and Kaya (2006). These parameters are shown 
below Table |4] Note that \/2k/3 = 0.6 < v = 1, hence the volatility process eventually hits 
zero with probability one, and the way the discretization schemes handle the zero boundary is 
strongly put to the test. 



As in Lord et al. (2010), we work with logS* instead of S to simplify matters. In other words, 
we consider equation (|3j) with d = 2, m = 1, and 



b{x) 



k{(3 



- xi, 

2 



a{x) 



for X = {xi,X2)~^ S E, and with X{0) = (wo,logso)^. Define (t{x) as the 2x2 diagonal matrix 
it(x) = diag[z^-^/Si, -^/xi], and S exactly as in (30). Then Condition Hi is satisfied if no > k and 
/ii satisfies < /ii < ^ {i^l^i'^ — ('^/^o)))^^^ (compare with Examplel? dealing with the CIR in 
isolation), and if the random vector e = (£1,62)^ is such that ei is a nonnegative variable with 
mean pi and variance 1, £2 is any variable with variance 1, and £1 and £2 have correlation p. In 
order to produce such a random vector for the simulations, we independently generate ei and £3 
as the Bernoulli-type variables with respective means pi = 0.657 and p^ = 1 and variance 1 (as 
in (20)), then we set £2= p ei + y^l — p^ £3. The payoff functional g{x) = e~^'^ {e^'^^'^^ —K)~^, for 
X = {xi,X2)~^ G C£;(M+), is continuous, as the payoff function of the corresponding put option 
would be. Convergence of the Monte Carlo prices to the right price when using our scheme 



in m«) hence follows from Proposition 10 and the put-call parity. 



All results in Table |4] are based on one million joint trajectories of the squared-volatility 
(following the CIR dynamics) and the log-price (i.e., N = 1000 000). No matter the scheme, 
the discretization of the log-price involes the square-root of the discretized CIR process, which 
might not be well defined if the discretized CIR becomes negative. The latter might happen if 



we use (bl), (b2) or (b4) on the CIR process. Following Lord et al. (2010), we fix this problem 
by using, for the discretization of the log-price, the positive part of the discretized CIR for 
schemes (bl) and (b2), and the absolute value of the discretized CIR for scheme (b4). Also, note 
that Kahl et al.l (2008) establish that the implicit Milstein method produces a well-defined and 



nonnegative approximate path for the CIR process only when the parameters satisfy ^/2k(3 > v. 
Our parameter values do not satisfy this condition, which is our main reason for not including 
such implicit schemes in this section. Indeed, with our parameters, the implicit scheme would 
have to be combined with another fix to work. However, the interested reader may consult 
Table 4 in [Lord et al. (2010), which is similar to Table |4] and includes results for the Milstein 
implicit scheme fixed to remain nonnegative. 

Table [4] confirms what has been observed in the bond pricing experiment: the first three 
schemes outperform schemes (b3) and (b4) in terms of biases. The first three schemes are also 
compared graphically in Figure [3] Lord et ah] (2010) numerically compare several discretization 
schemes for option pricing in Heston's model and conclude that scheme (b2) is very efficient. 
From Table |4] and Figure |3j one can clearly conclude that our scheme is competitive and can be 
compared favourably with scheme (b2). 
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Table 4. Pricing of a European call option 



/fc — o Lcpt) / y cdl 


I— < /Z\T* 11 1 1 1 1 1 

Ocl llOUlil 










5 


-0.121 

I^U. iUo j 


1.868 


0.359 


8.318 
/'n 1 Q/i ^ 


6.995 

I^U.ioo j 


10 


— U.Uo / 

(0.111) 


n OAS 
(0.120) 


U. ioO 

(0.115) 


D.UOO 

(0.165) 


41.40O 

(0.158) 


20 


-0.061 

(0.0.112) 


0.500 
(0.116) 


0.137 
(0.113) 


4.419 
(0.148) 


2.733 
(0.140) 


40 


-0.070* 

(0.112) 


0.147 
(0.115) 


-0.013 

(0.114) 


3.181 
(0.139) 


1.649 
(0.129) 


80 


-0.013* 

(0.113) 


0.126 
(0.114) 


0.057 

(0.113) 


2.377 
(0.131) 


1.074 
(0.123) 


160 


0.093* 

(0.114) 


0.059 

(0.114) 


0.030* 

(0.113) 


1.795 
(0.127) 


0.651 
(0.119) 


rate of conv. 


0.50 


1.01 


0.92 


0.45 


0.69 



Bias and margin of error at the 95% confidence level (in parentheses). 
Volatility parameters: k = 2, (3 = vq = 0.09, u = 1. 

Option and stock price parameters: T = 5 years, sq = 100, K = 100, r = 0.05. 

Correlation parameter: p = —0.3. 

True option price: 34.9998. 

Means for £i and £3: /^i = 0.657 and //3 = 1. 



Figure 3. Bias comparison for the option price 
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As in Section 6.1 orders of convergence are estimated. Rough estimates in Table |4] may be 
compared with those in Table 4 in Lord et al. (2010), who claim that it is quite hard in this case 
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to properly estimate the order of convergence even when using 10 milhons trajectories. However, 
care is advisable when interpreting the estimates of the order of convergence in Table [4j especially 
for our method. Indeed, the biases presented in Table [4] are not significantly different from zero 
(for our scheme only) for values of n as small as 10 and 20; moreover, the bias for n = 5 is the 
lowest. A more accurate estimate is therefore provided in Table [5] Although the estimated rate 
of convergence of 0.65 is smaller than those estimated for methods (bl) and (b2), it has to be 
kept in mind that (as is made clear in Figure |3]) the biases for the proposed method are small 
overall in comparison with other methods. Hence, one may certainly claim that the proposed 
method is a competitive alternative to other methods, even more since Figure [T] and Figure [2] 
led to a similar conclusion in another meaningful experiment. 



Table 5. Option pricing using the Bernoulli distribution 



n 


N 


bias 


margin 


RMSE 


5 


5 X 10'' 


-0.1144 


0.0480 


0.1169 


10 


10 X 10^ 


-0.0911 


0.0350 


0.0929 


20 


20 X 10^ 


-0.0435 


0.0251 


0.0453 


40 


40 X 10*^ 


-0.0452 


0.0178 


0.0461 


80 


80 X 10^ 


-0.0227 


0.0126 


0.0236 


160 


160 X 10^ 


-0.0109 


0.0090 


0.0118 


rate of conv. 


0.654 




0.641 



Bias, margin of error at the 95% confidence level, and RMSE 
for different combinations of n (the number of steps per year) 
and (the number of trajectories). 
The parameters are the same as in Table [4] 



Finally, Tables [T] [2] and [4] clearly illustrate that the rate of convergence, as measured here, is 
not as good a measure of precision as it appears to be. In these tables, the error of the biases 
are comparable but our proposed method is much more precise in the sense that the bias is not 
significantly different from zero, even for small number of steps n. 



6.3. Conclusion. In summary, in addition to the simplicity of our scheme, its great flexibility 
and the fact that approximate prices are theoretically known to converge, it is observed to be 
numerically competitive with other existing schemes in two representative experiments. We 
recall that our scheme and our convergence results hold in more general diffusion models and for 
a wide variety of payoffs. 
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Appendix A. Proof of Theorem El 



We now prove Theorem [S] As discussed in Section [3j we will use the martingale problem 
formulation. We need the next result, which gives sufficient conditions for Anf defined in ( 15 ) to 
converges to Af, for all / E C^(M'^), and then concludes that this is enough for the corresponding 
processes to converge as well. Throughout, | • | is the Euclidean norm. 

Proposition 16. Assume Conditions^ and^ For each n > uq, let Kn{t,x,r) be a time- 
dependent transition function defined on x E x B(R'^) and such that Kn{t,x, E) = 1 for all 
{t, x) G M+ X E. For n > no, set 



(31) 

and 
(32) 



hn{t,x) := n {y - x)Knit,x,dy) 

J\y-x\<l 



Q,fi (tj x) .- 



n\ {y - x){y - x)'^ Kn{t,x,dy) , 

'\y-x\<l 



for each (t, x) G M+ x E. Assume further that, for any r > and e > 0, the following sequences 
tend to zero as n goes to infinity: 



(33) 

and 
(34) 



sup \bn{t,x) 
\it,x)\<r 



b{t,x)\, sup \an{t,x) — a{t,x)\ , 

(t,x)eR+xE 
\it,x)\<r 



sup nKnit,x, {y : \y - x\ > e}) . 

{t,x)£R+xE 
\it,x)\<r 



Let {Yn{k))k>o be a Markov chain, with 1^(0) = xo G E, and transitions governed by Kn through 
equation (14). Then 

lim sup \Anf{t,x)-Af{t,x)\=0, 

(t,x)m+xE 

\{t,x)\<r 

for all r > and f G C^(M°'), where An is defined in terms of Kn in (|15|) and A is defined 



in (11). Moreover, the sequence of processes (^n)n>no; defined from (^^)n>no 

through 

converges in distribution to the solution of the martingale problem for {A,xq). 



That proposition corresponds to results in (Stroock and Varadhan [1979 Section 11.2), re- 



stricted to E. See also Corollary 7.4.2, in conjunction with Proposition 3.10.4 in Ethier and 



Kurtz (1986). Note that the functions a and b postulated in Section 11.2 in Stroock and Varad 



han (1979) or Corollary 7.4.2 in Ethier and Kurtz (1986) are time-homogeneous. To extend these 



results to the case where a and b are time-dependent, one only needs to append the time to the 
state variable. In other words, consider the ((i-l-l)-dimensional state process X{t) := {X{t)'^ ,t)'^ . 
Then X is a solution to the SDE Js]) if and only if X is a solution to dXt = b{Xt)dt + a{Xt)dW{t) 
with initial condition X[0) = (xq70)^, upon defining 



where W stands for a {d+ l)-dimensional standard Brownian motion. 



a{t, x) 
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Our goal is now to verify that Proposition [16] applies when the transition function Kn is the 
one defined at ( 13 ). From Condition |4| it is clear that Kn{t, x,E) = 1 for all {t, x) G E. The two 
following lemmas establish that the three sequences given in (|33|) and ( 34 ) tend to zero as n goes 



to infinity, which is all that is needed to get the desired convergence. To lighten the notation, 
we shall write Z := e — fj,. 

Lemma 17. Assume Conditions^ andj^ For t > 0, x G E and e > 0, define 



Dn{t,x,e) :-- 



—b(t, x) H — i=^(t, x)Z 
n Jn 



> e 



Then, for any r > and e > 0, we have, as n goes to infinity, that 



(35) 
(36) 

Proof. For each r > 0, define 
M,(r) 



sup E [\Z\ lD„{t,x,e)] 
{t,x)eR+xE 
\{t,x)\<r 

sup nF{Dn{t, X, e)) 

{t,x)£R+xE 
\{t,x)\<r 



0: 



0. 



max |6(t,x)| and M^{r) 
{t,x)eR+xE 
\{t,x)\<r 



max |(j(t, x)| . 

{t,x)eR+xE 
\(t,x)\<r 



Since b and a are assumed to be continuous, we have < Mf,(r), M^(r) < oo for every r > 0. Fix 
r > and e > 0. For {t,x) eR+ x E satisfying \{t,x)\ < r, note that \^b{t,x) + -^a{t,x)Z\ < 

^Mi,{r) + ■^M^{r)\Z\, from which one may establish existence of 7 > and Nq > (depending 

on r and e) such that 



(37) Dn{t,x,e) C {\Z\ > ^/nJ} , for all n > iVo and |(t, x)] < r . 

If |(t, 2;)| < r and n > Nq, we then have 



{|Z|>v^7} 



E[\Z\HD„it,x,e)] <IE 

and since the right-hand side does not depend on (t, x), ( |35[ ) easily follows from the dominated 
convergence theorem, and the fact that the components of Z have finite second moments. The 



inclusion in (37) also yields the first inequality in 

I 



n 



1 



'n{t,x,e)) < -^{Vn-fYE 



T 



{|Z|>v^7} 



< 



1 



{|Z|>v^7} 



for |(t, x)| < r and n > Nq, and (36) is also immediate from the dominated convergence theorem. 



Note that (36) is equivalent to saying that, with Kn defined in (13), the supremum in equa- 



tion (34) tends to zero as n goes to infinity, for all e > and r > 0. 



Lemma 18. Assume Conditions and^ For the transition function Kn defined in (13), 
and On and bn defined in (32) and (31) respectively, the two sequences in (33) tend to zero as n 



goes to infinity, for all r > 0. 
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Proof. From (31) and E[Z] = (recall Condition |4j), we can write 

"1 



bn{t,x) 



nE 



b{t,x) + -^a{t,x)Z ] I 



= b{t, x) - b{t, x)¥{Dn{t, X, 1)) - V^a{t, x)E • 

Here, D„ is as in Lemma [17] Using the Cauchy-Schwarz inequality, and with Mh and M5- as in 
the proof of Lemma 17 we then see that, for < r. 



\bnit,x)-b{t,x)\ < Mbir)F{Dn{t,x,l))+M^{r) (E\Z 



2\ 2 



nF{Dn{t,x,l)) 



The result follows by ( 36 ) in Lemma JT7| Similarly, note that 

1 



an{t,x) 



nE 



—b(t, x) H — i=o'(t, x)Z 
n \ n 



Rearranging, and using the fact that E\ZZ'^\ 
an{t, x) — a{t, x) 



S and a = ctScj^ (see Condition U]) , we get 



n 



-b{t, x)b{t, x) 



,{t, X, 1)) - a{t, x)E[ZZ ' Ii,„(t,,,i)]a(t, x) 



+ ^6(t,x)]E[zTl^^-^]a(t,x)T + ^a{t,x)E[Zl^^^^]b{t,x) 



If x)\ < r, we then have 

\an{t,x)-a{t,x)\ < -Mb(r)2 + AM5(r)Ms(r)E|Z| +Ms(r)2E[|Z|%„(t,^,i)], 



and the result follows from (35) in Lemma 17 
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